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1. Introduction 

Driving a system away from thermal equilibrium requires work. Since thermal 
fluctuations play a paramount role in small, mesoscopic systems it is only consequent 
to also define a fluctuating work along single trajectories [TJ[2]. While the celebrated 
Jarzynksi [3| and Crooks |4; i5j non-equilibrium work relations constrain the possible 
shapes, the actual distribution of the work is still a non-universal function that 
depends on both the system dynamics and the driving protocol. Analytical expressions 
for the work distribution in isothermal processes are rather scarce except when the 
distribution is exactly a Gaussian [B]. In particular, a single Brownian particle in a 
moving harmonic potential has been studied extensively both theoretically [71 |H1 El H] 
and experimentally [101 lllj by trapping a colloidal particle with optical tweezers. 

The direct application of the Jarzynski relation to numerical or experimental 
data in order to extract equilibrium free energy differences is marred by the fact 
that it falls into the class of biased estimators. Rare, untypical trajectories with 
low work values are exponentially weighted and a large number of observations is 
required for estimates to converge. This number typically grows exponentially with 
the mean dissipated work. Two principal schemes have been discussed to address 
this problem: (i) reduction of the mean dissipation by using optimal protocols [12] 
or escorted simulations (13) : (ii) unbiased estimators based on Bennett's acceptance 
ratio method |14j or extended bridge sampling [T5] . 

Of course, the easiest method to reduce the mean dissipated work might still 
be to reduce the driving speed and keep the system close to equilibrium. In this 
case a Gaussian has been predicted [El [IT] as the limit distribution for the work. 
However, as has been noted many times, a naive application of a truncated Gaussian 
work distribution leads to wrong results for the free energy difference (see, e.g., the 
discussion in Ref. ^18j). While the central limit theorem predicts that the center 
of the distribution approaches a Gaussian this does not necessarily apply to the 
extreme tails of the distribution, which are crucial for the correct determination of 
free energy changes. For practical purposes it is, therefore, important to understand 
the convergence of the work distribution towards its limiting Gaussian distribution. 

In this paper we study the probably simplest system that leads to a non- 
Gaussian distribution of work: a single Brownian particle moving in a one-dimensional, 
tightening harmonic potential. In addition to an exponential tail for large work values 
the probability for negative work values is exactly zero. Instead of determining the 
distribution of work directly we consider its moment generating function and obtain a 
closed set of non-linear ordinary differential equations. Using this exact result allows 
us to study the work distribution as we reduce the driving speed. 

2. Work distribution and its generating function 

We consider a particle with position x moving in a harmonic potential 



with time-dependent strength k{t). Assuming overdamped dynamics the evolution of 
the probability distribution pQ{x,t) is well described by the Fokker-Planck equation 




(1) 



dtpo = Ht)po, 



L{t) = d^[k{t)x + d,]. 



(2) 
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Here and throughout the paper we measure energy in units of k^T with Boltzmann's 
constant k^, length in units of ^/k^^Tjk^, and time in units of k^T /[D^ko). The free 
diffusion coefficient of the particle is , and ko is the initial non-dimensionless value 
of the trap strength (i.e., from now on fc(0) = 1). Given a driving protocol k{t) from 
t — to t — ts the work is defined as the functional 

w[x{t)] = dt -^{m,t) ^-j^ dt k{t)x\t) (3) 

with switching time t^. In particular, for /c > (fc < 0) the work is always non-negative 
(non-positive). 

Since x is the result of a stochastic process the work also will be a random variable 
with distribution p{w). For any driving protocol k{t) the Jarzynski relation 

(e-"-) = y" d«; e^Xu;) = er^'" (4) 

provides the link between non-equilibrium work and the equilibrium change of free 
energy Ai^, which for the system studied here becomes 

AF = F(i,)-F(0) = ilnfc(i,). (5) 

For the analytical study of work distributions it has turned out to be convenient to 
consider the joint probability p{x, w, t) to find the system in state x at time t and to 
have spent the accumulated work w so far. Its time evolution is governed by 

dp dU dp 1 ; 2 5^ 

d-t^''P-^d^^''P~2'''' d^- 
In the following we work with the Laplace transform 

nOO 

px{x,t)= / dw e'^'"p{x,w,t). (7) 

JQ- 

For simplicity we focus on fc > 0. After one integration by parts and using 
p{x, 0^ ,t) ~ we obtain from ([6| the evolution equation 

dtpx Lpx - ^x'^px- (8) 

This is a linear reaction-diffusion type equation, called a 'sink' equation, where 
diffusion is governed by the operator L and 'probability' is created or annihilated 
with a rate proportional to A, i.e., the function px is not normalized. Rather, its 
integral is the moment generating function 

Mt) = {e-^n = J dxpx{x,t). (9) 

In particular, the mean and variance are obtained as 

{w)^~pi, al = {up) - {wf ^ P2 - pi, (10) 

where pn = d"?/'A/dA"|A=o- 

Already in one of the first experiments demonstrating a non-equilibrium work 
relation |19| it has been noted that the distribution p{w) for the work is distinctly 
non-Gaussian. More recent computer simulations have found exponential tails for 
large w [20j [21] . Indeed, considering the extreme case of instantaneously switching 
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from fc(0) = 1 to fc(ts) = 1 + A at some arbitrary time < t < tg the work spent is 
Ax'^(t)/2. The work distribution, therefore, becomes 

where we have averaged over the equihbrium distribution 

Peq(a;) - (27r)-l/2e-5-^ 
The generating function reads 

ip^ = dwe-^Xw) = (l + AA)-i/^ (12) 



0- 



.-1 



As a signature of the exponential tail -0^° diverges for A — > A* with A* = —A" 

The other limiting case is that of a quasi-static process with — > cxd and work 
distribution p'^^{'w) = S{w — AF), implying the generating function 

V^f (t) = e-^'^^ = = [Ht)]'^^^- (13) 

As a third exact result we note that for A = 1 

p,{x,t) = (27r)-V2e-^M*)-=, ^^(t) ^ e-AF (^4) 

solves the sink equation ([8|. The generating function tpi obeys the Jarzynski relation, 
cf. ^ with Q. This particular solution to the sink equation Q has been discussed 
first by Hummer and Szabo [22j . 

3. Time evolution of the generating function 

From Eqs. ^ and ^ we obtain the equation of motion 

after integration over x, where 

cj}x{t) = J dxx^pxix,t) (15) 

is the generalized second moment. Multiplying ^ by and integrating again over 
X results in 

Xk f 

4>x = -2k<f)x + 2ipx ^ ^ / da; x'^px- 

By following this scheme we obtain a hierarchy of coupled ordinary differential 
equations. Fortunately, for the harmonic oscillator a closure can be found as follows. 



For A = 1 we observe that the solution ( 14 ) is a, albeit not normalized, Gaussian 



Since otherwise A = 1 is not special we conclude that px is a Gaussian for all A. Then 
the closure 

„4 ,, /"^^ ^4 PA _ 3(/>2 



dxx^A-V^A / dxx^^^^ (16) 

follows from Gaussian statistics, expressing the generalized fourth moment in terms 
of the functions tf^x and (j)x- The resulting system of non-linear, first-order ordinary 
differential equations 
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Figure 1. The generating function VaC* = ^s) A for A = 2 and different 
driving speeds re. Tfie shaded regions indicate the accessible range bounded by 
the quasi-static | |13| l and instantaneous ( |12| hmiting cases (dashed hues). All 
lines cross at A = (normalization) and A = 1 [Jarzynski relation The right 
panel shows the decay of ipx for large A. 



is our first main result. Through inserting the Gaussian ansatz 



■ exp 



x^^x{t) 



^ 2T:(bx{t) t 2(j)x{t) 

into equation ([s]) it is straightforward to check that (17) is indeed correct. Augmented 
by the initial conditions tp\{0) = <j>\{0) = 1 these equations are readily solved 
numerically by standard techniques. 

In figure [l] we plot the generating function i^x{t — tg) at the final value of the 
control parameter for different driving speeds k = k using the linear protocol 

k{t) = l + Kt, Kts = A. (18) 

The second relation fixes the switching time t^. In the right panel of figure [T] the 
generating function Tp\ for large A is shown. The probability p(O^) to have spent no 
work is related to this asymptotic behavior through the initial value theorem, 

p{0+) = lim AV-A- 

A— ^oo 

The numerical results suggest that even for large driving speeds Va decays faster 
than A~-^ with ^(0"*") = and only for instantaneous switching p{0~^) oo, see 
equation (111. For large w an exponential tail is expected [H], which corresponds to 
a divergence of ipx at A* < 0. Starting from ipxi^) = 1 we expect a singularity at time 
t* for all A < A*, where A* corresponds to t* = ts- Together with the results Eqs. ( [T2| ) 
and ( 13 ) for the two limiting cases this implies a value for A* that moves from — A^^ 
to —00 with decreasing driving speed k. 



4. Slow driving 

Of greater practical importance than quasi-static driving with k — > are processes 
with slow but finite driving speed k. In this case a Gaussian distribution for the work 
has been predicted [T7]. Of course, in the present situation the work distribution can 
never be strictly a Gaussian since ^(w ^ 0) = and because of the existence of an 
exponential tail for large w. 
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Figure 2. a) Difference between true free energy AF and the truncated 
approximation ( |19[ l for the linear protocol ( |18[ |. The dashed line shows 
convergence with ~ k^. Both mean work b) and variance c) grow linearly with 
driving speed (dashed lines ~ k). 



We calculate the true mean and variance [see ( 10 1] by integrating the four coupled 
equations 



Ail 



Ai2 



= -2k(j)n + 2, 



(1) - 



'2k(j) 



(1) 



2^1 - 



3k 



which follow from ( 17) after expanding tpx ~ I + /iiA + /i2A^/2 and (px ~ (po + (/''■^■'A. 
Truncation of the cumulant expansion for the free energy leads to the approximation 

1 



AFe ^ (w) - -ai 



(19) 



In figure [2ji) the difference AF — AFc is plotted as a function of the driving speed 
n for the linear protocol (18). It shows that the difference vanishes as and that, 
therefore, the true free energy can be arbitrarily well approximated by equation ( 19 ) 
through lowering the driving speed k. In figure [2]d) and c) the mean dissipated work 
and variance are shown, respectively. 

The strategy used in Ref. |17] to obtain the work distribution for slow driving 
corresponds to expanding the generalized second moment 



(20) 



in powers of k. To first order we obtain from (17 1 i\)x ~ ~AK'i/'A/(2A:), which is solved 
by the quasi-static solution ilf^ (13). For the next order we plug the expansion (20) 

into ( |l7| and retain only terms of order k with (^^^^ ^ k. The result is 



-(l + V2)|^-2fc0«-§^. 



with 



1 - A 



The generating function now reads 
Wx\ts) = exp 



Infc 



A(l-A), 



(21) 
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Figure 3. Comparison of the initial slope of mean [figure [2J3)] and variance 
[figure [2|;)] with equation l |22| . The solid line shows 7 and the dashed line shows 
7/2. 



The corresponding distribution is of course a Gaussian with variance — 7K and 
mean (w) = AF + 7^/2, where 

In contrast to the truncated Gaussian distribution obtained by using the true mean and 
variance, the distribution foUowing from (21 1 always fulfills the Jarzynski relation Q. 
In figure |3] we compare the initial linear slope of mean and variance for these two 
Gaussian distributions, where the slope of the numerical data has been determined 
through a polynomial fit. The agreement between the two Gaussian distributions is 
our second main result. 



5. Conclusions 

For a driven harmonic oscillator we have derived the equation of motion for the 
moment generating function of the work distribution. Even though the exact 
distribution is non-Gaussian it can be approximated by a Gaussian such that the 
error vanishes with the square of the driving speed k^. While such a behavior is 
expected to hold in general [T^ it remains to be investigated to which extent the 
exponent 2 is system dependent. We have focused on a linear protocol with fc > 
which corresponds to a stiffening trap. The case A; < follows in analogy by defining 
the Laplace transform in ([7| with respect to negative work values. However, it should 
be kept in mind that slow driving is defined with respect to the time-scale separation 
between driving speed and the system's relaxation time 1/k, which for a widening 
trap increases strongly. 

The extension to a time-dependent quadratic potential energy with many degrees 



of freedom is straightforward through using the Wick theorem in (16). It will 
be worthwhile to study the method employed in this paper for more complicated 
potentials U{x). While in the present case the Gaussian closure is exact, approximate 
closures for other potentials might nevertheless lead to accurate results for the work 
probability and cumulants. 
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